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Abstract 

A variety of problems in acoustic and electromagnetic scattering require the evaluation 
of impedance or layered media Green's functions. Given a point source located in an un¬ 
bounded half-space or an infinitely extended layer, Sommerfeld and others showed that 
Fourier analysis combined with contour integration provides a systematic and broadly ef¬ 
fective approach, leading to what is generally referred to as the Sommerfeld integral repre¬ 
sentation. When either the source or target is at some distance from an infinite boundary, 
the number of degrees of freedom needed to resolve the scattering response is very mod¬ 
est. When both are near an interface, however, the Sommerfeld integral involves a very 
large range of integration and its direct application becomes unwieldy. Historically, three 
schemes have been employed to overcome this difficulty: the method of images, contour 
deformation, and asymptotic methods of various kinds. None of these methods make use 
of classical layer potentials in physical space, despite their advantages in terms of adap¬ 
tive resolution and high-order accuracy The reason for this is simple: layer potentials are 
impractical in layered media or half-space geometries since they require the discretization 
of an infinite boundary In this paper, we propose a hybrid method which combines layer 
potentials (physical-space) on a finite portion of the interface together with a Sommerfeld- 
type (Fourier) correction. We prove that our method is efficient and rapidly convergent for 
arbitrarily located sources and targets, and show that the scheme is particularly effective 
when solving scattering problems for objects which are close to the half-space boundary or 
even embedded across a layered media interface. 

1 Introduction 

Problems of acoustic and electromagnetic wave scattering in half-space or layered media ge¬ 
ometries require the solution of the the governing partial differential equations subject to suit¬ 
able boundary and radiation conditions. In the two-dimensional, time-harmonic setting, both 
reduce to the Helmholtz equation (Figure 

Au + k^u = /, (1.1) 
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Figure 1: (a) Scattering in the presence of an impedance half-space, with a point source 
defining the incoming field and a sound soft scatterer Oq. In (b), the impedance boundary is 
replaced by two-layer media, with a distinct Helmholtz parameter in the lower half-space. 
The scatter Oq is partially buried and touches both media. 


with boundary conditions enforced on a scatterer Oq and interface conditions enforced on the 
line y — 0, either of impedance (Robin) type 


X — h ikocu = 0, 
an 


( 1 . 2 ) 


or of transmission type 


[u] = 0, 


du' 

dn 


(1.3) 


The Helmholtz coefficient k is given ask — w / c, where m is the governing frequency and c is 
the wave speed in the medium. For the sake of simplicity, we will assume that on the scatterer 
Oo the total field u satisfies homogeneous Dirichlet boundary conditions 


^ — 0|ano- 


In electromagnetics, this condition corresponds to the case of scattering from a perfectly 
conducting obstacle in transverse-magnetic (TM) polarization, and in acoustics to the case of a 
sound-soft obstacle. Here and in what follows, n = (0,1) is the unit normal on the line y = 0, 
denotes the partial derivative of u in the normal direction, and a is an impedance constant fP] 
with 5ft(a) > 0. The expression [/] denotes the jump in the function / across the line y — 0, 
which we will denote by T in the remainder of the paper. We will denote by the incoming 
field induced by the sources / in 0- We will limit our attention, without loss of generality, 
to either point sources or plane waves. To ensure uniqueness of the boundary value problem, 
a radiation condition must be imposed to enforce that the scattered field is decaying. Thus, 
we assume that the total field u is written in the form u — u' + w®, where the scattered field w® 
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satisfies the Sommerfeld radiation condition |TT| : 


lim Vr (^ - iku^^ 

r^oo \ dr J 


0. 


(1.4) 


We will also assume that the Helmholtz parameter k is constant in either the upper or lower 
half-space, with 3ft(A:) > 0 and ^(k) > 0. Some applications require variants of the interface 
conditions above, such as 



(1.5) 


where 7,^6 are piecewise constant material parameters |[T0ll22ll24ll . The method of this paper 
extends to these cases in a straightforward manner. 

Integral equations are natural candidates for solving the problems described above since 
they discretize the scatterer alone and impose the Sommerfeld radiation condition by con¬ 
struction. In order to make effective use of this approach, however, one must generally evalu¬ 
ate the governing Green's function which satisfies the homogeneous interface conditions ( | 1 . 2 [ ) 
or ( |1.3[ ). This avoids the need to discretize the interface T (the infinite line y = 0), and the 
most common treatment, using Fourier analysis, was pioneered by Sommerfeld, Weyl and 
Van der Pol ||34l|3ZHMI- Over the past several decades, a number of methods have been de¬ 
veloped, based on this approach. These include using ideas from high-frequency asymptotics, 
rational approximation, contour deformation (ZlElEZHsTJI^, complex images ||26[|35H36l, and 
methods based on special functions or physical images Il 2 ^ . 

Without reviewing the methods listed above in detail, we simply note that all of them are 
aimed at the efficient pointwise evaluation of the impedance or layered media Green's function, 
rather than treating the scattering problem itself in a more unified fashion. Here, we consider 
a substantially different approach, motivated by the fact that close-to-touching interactions be¬ 
tween compactly supported scatterers are easily accounted for using standard physical space 
layer potentials, which can be adaptively refined to cope with the near singularities induced by 
the geometry (see, for example, ||l3[[T6lllZ|). Potential theory cannot be used naively in layered 
media, however, because of the infinite extent of the interface. In principle, however, it seems 
plausible that the failure of rapid convergence of the Sommerfeld integral is due entirely to the 
close-to-touching interaction. This turns out to be the case, and we present a rigorous, hybrid 
method for scattering problems that combines the best features of layer potentials (adaptiv¬ 
ity with high order convergence) and of the Sommerfeld representation (spectral accuracy for 
smooth functions). Put differently, layer potentials in physical space will be used to capture 
features of the solution with high-frequency components in the Fourier domain, and the Som¬ 
merfeld integral will be responsible only for the remaining low-frequency components. For 
this we solve a local integral equation and apply a smooth window function to the solution in 
order to capture the most singular behavior. 
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Remark 1. It is worth emphasizing that window functions have been used previously to ac¬ 
celerate the evaluation of layered media Green's functions IZl IH |3TJ |32l. The approach in 
those papers, however, is based on using carefully chosen partitions of unity and asymptotic 
analysis to evaluate slowly decaying and oscillatory Sommerfeld integrals for each source and 
target. We are using a partition of unity in physical space to enforce rapid decay in a single 
Sommerfeld-type correction that can be used for all target points. 

Bruno and collaborators ISJ |6| have independently developed a method that makes use of 
the same underlying intuition. They also propose using a window function regularizing the 
integral operator in physical space in order to handle the most complicated features of the 
scattering problem. Unlike our scheme, asymptotic methods and stationary phase arguments 
are used to approximate the solution of the integral equation on the entire interface. Their 
windowing scheme could be adapted for use in place of our local integral equation (described 
below), with the potential for eliminating the artificial endpoint singularities introduced in our 
scheme. This would reduce the complexity of the implementation |^. 

An outline of the paper is as follows: in Section]^ we review some basic properties of the 
free-space Green's function, layer potentials, and their spectral representation. In Sections 
and 1^ we focus on the construction of a local integral equation and describe our hybrid rep¬ 
resentation in more detail. We then prove, in Section that our local integral equation is 
well-posed. In Section]^ we prove decay estimates on the integrand of our Sommerfeld correc¬ 
tion. Sectionj^contains some illustrative examples for both pointwise evaluation of the Green's 
function and for scattering computations in the presence of obstacles. Section [^contains some 
brief, concluding remarks. 


2 Spectral representation of the Green's function 


For A: G C with non-negative imaginary part, it is well-known that the Green's function for the 
free-space Helmholtz equation ( |1.1[ ) is the zeroth order Hankel function of the first kind: 

Gkix,xo) = - acol). (2.1) 

It satisfies 

{A + k^)Gk{x,xo) = S{x - xo), ( 2 . 2 ) 

where ^(x — xq) represents the delta function at xq. The Somm erfe ld integral representation 
for is obtained by taking the Fourier transform of equation \2.2\ with Fourier coordinates 
(Ax/ Ay), and evaluating the integral in Ay via contour integration This yields 



g-VA2-F|y-yo| 

VA^-k^ 


^iX(x-xo) 


(2.3) 
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Since Aj is the only Fourier variable remaining after contour integration, we denote it simply 
by A. 

The formula ( |2.3[ ) plays a central role in scattering theory, and numerous schemes are avail¬ 
able that yield high order accuracy upon discretization. The simplest, perhaps, is contour de¬ 
formation to avoid the square root singularity in the integrand (see for example, EKTOl). We do 
not review the literature here, since we are primarily concerned with the range of integration 
required in the Fourier integral parameter A. The relevant considerations are most easily under¬ 
stood in the context of computing the Green's function for a perfectly conducting or sound-soft 
half-space, where we seek to impose the Dirichlet condition u — 0 on the interface F. 

The obvious solution, of course, is to construct the corresponding Green's function, which 
we denote by G^, using the method of images. Assuming xq = (xo,i/o) with yo > 0, let its 
reflected image Xq — (xq/ “J/o)* h is then easy to verify that 


G^{x,xo) = ^-H^^\k\x - xo\) - ^-H^^\k\x - x^\). 


(2.4) 


Using the Sommerfeld integral, we may instead write 


Gt{x^xo)='^Hi^\k\x-xo\)-^l_ 


00 p-V^^-k'^\y-\-yo\ 


VA2 - 


^/A(x-xo) 


(2.5) 


While this is more complicated than ( |2.4[ ), the analogous approach can be used even when the 
boundary or interface condition does not support a simple image representation, as we shall 
see below. For the moment, we simply wish to observe that when |y + yo | is large, the integrand 


in (2.5) is rapidly decaying, once 3ft(A^ — k^) >0. When |y + yo| is small, however, the decay is 
negligible and the range required in A is large. In this regime, the Sommerfeld representation is 
extremely inefficient unless various asymptotic or contour deformation methods are employed. 

A third approach to imposing the homogeneous Dirichlet interface condition u — 0 on T 
is to write u{x) = {k\x — xo|) + u^{x), with the scattered field represented as a double 
layer potential with unknown density a: 


u^x) = D’^[(t]{x) = ^ 


dn 


jGkix^x') 


cr(x') ds{x') 


( 2 . 6 ) 


Here s denotes arclength along T, and d/dn' denotes differentiation in the outward normal 
direction at the point x'. Substituting the representation of w® into the boundary condition 
M = 0 yields a second-kind integral equation; 

^ + D’^*[cr] = -u\ (2.7) 

where Dp* denotes the principal value of Dp. On a half-space, however, it is well-known that 
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Figure 2: (a) For a point source at (0,2) with k = 1, we plot the imaginary part of the 
function on F. The response is smooth, oscillates over the indicated range (—100,100), 
and decays slowly (like 1/ y/x). (b) When the point source is near F, the oscillations and 
slow decay are still present, but the function is nearly singular at the close-to-touching 
point. 


Dp* vanishes IITSlI , so that a — —2u^ and 


u^(x) = 


-2u' 


h)- 


This is not useful in practice, because T = (— 00 , 00 ) and w' is slowly decaying in space, so that 
the range of integration is extremely large (see Figure]^. Note, however, that when yo is small, 
w' is nearly singular only at the close to touching point in physical space. This accounts for the 
slow convergence of the integrand in the Sommerfeld integral, which is simply 

the Fourier transform of u' along the line F. 

Our hybrid approach is based on the following premise; that the poor convergence of 
the Sommerfeld integral is due entirely to the near singularity in the layer potential density 
{cr = at the close-to-touching point. Thus, we partition w' into a local part, for which we 

may make effective use of layer potentials, and a remainder, for which we can effectively use 
the Sommerfeld representation (see Figure]^. 

More precisely, let us assume we have a window function W(x) 6 C~(]R), supported on a 
finite section of the interface, Fq = (—Mq, Mq), that satisfies 


lV(x) 


0 forx<—3/4Mo, 

1 for — l/4Mo < V < l/4Mo, 
0 for X > 3/4Mo. 


( 2 . 8 ) 


We can construct such a window function from a compactly supported C°° function, such as 
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Figure 3: For a point source near the boundary, the double layer density 2u^ is nearly sin¬ 
gular on F. We use a smooth window function (see text) to capture the nearly singular 
part in physical space (top right), and the Sommerfeld integral to account for the remainder 
(bottom right). 


the standard hump function: 


Y(x) 


e for |x| < 1; 

0 otherwise. 


For this, we let <I>(x) = ^ fj Y(t) dt, where A — JgY(t) dt, and define 

W(x) = Wmo( x) = l(<I)(x + Mo/2)Mo/2)) . (2.9) 

It is straightforward to verify that, once Mo > 4, W(x) satisfies the desired conditions. We use 
notations W and Wmo interchangeably, depending on the particular relevance of the parame¬ 
ter Mq. 


Remark 2. In practice, we use 

W(x) = Wmo(x) = 1/2 (erf(x + Mo/2) - erf(x - Mo/2)) (2.10) 

where erf(x) is the error function 


erf(x) = dt. (2.11) 

y/n Jo 

When Mo = 28, W in (|2.10|| satisfies (|2.8|l with more than fourteen digits of accuracy. For 
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smaller windows (smaller values of Mq), W can simply be rescaled. The smaller the window, 
however, the less rapid the decay of its Fourier transform. 

Letting crw(x) = lV(x) • cr(x) = Wm(,(x) • (—2w'(x)), we seek to represent the scattered field 
in the form 


= / 
Jy 


r. 


c^w(x')ds(x') + L I” 


( 2 . 12 ) 


Using the decomposition a = Cw + (1 — W) (7, it is straightforward to verify that 


^w(A) = -2(l-W)* 




-zAxo 


(2.13) 


where /(A) denotes the Fourier transform of /(x): 

/ oo 

/(x) e-'^^dx. 

-OO 

That is, w® is represented as a double layer potential over a finite region To plus a Sommerfeld 


correction with density defined in (2.131. It turns out that is rapidly decaying as a 


function of A, as shown in Theorem 2.2 below. In subsequent sections, we show how this 
hybrid representation can be applied to the cases of interest, rather than merely the Dirichlet 
problem where one would (of course) use the simple image solution in practice. 


Lemma 2.1. The Fourier transform of the window function 'Wmq in ( |2.10 ) decays super algebraically. 
Proof This follows immediately from the fact that W is C°° and integration by parts. □ 

Theorem 2.2. Let w® be the scattered field induced by a point source at {0,h) satisfying the Dirichlet 


boundary condition w® = —uL If we represent w® in the form (2.12), then |fw('^)| decays superalge- 
braically, independent ofh. 


Proof. Rather than deriving this estimate directly from formula (2.13|l, we write 


lw(A) = —2 


1 — Wmo * w* 


= -2 


1 U^Mn + ^Mn) 


where denotes the field due to a smooth compactly-supported distribution centered at the 
source point (0, h) of the form 

^Wi3\x-i0,h)\). 

This scaled version of W vanishes at a distance Mo/4 from the center. If we let 

rMo/ 4 : 

B = 2n Jo{kp)W{3\x-{0,h)\)pdp, 

Jo 
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then and are identical, once \x — {0,h)\ > Mo/4. This is easily established by observing 
that for |x| > Mo/4, the multipole expansion induced by the smooth, compactly supported 
source and the point source are identical. The choice of B follows directly from the Graf addi¬ 
tion theorem Il28ll . Since 1 — Wmq is identically zero in the interval [—Mo/4,Mo/4], it is clear 
that the product (1 — Wmq) identically zero on the real axis. The result now follows 

from Lemma [ 24 ^ and the fact that (A) decays superalgebraically, independent of h. □ 


Definition 2.3. Let the functions /,g : R ^ C. The functions / and g are said to be spectrally 
equivalent in the far-field if (/ — g) (1 — W) =0. 


Thus, Theorem 


2.2 


relies on the fact that and u 


Mo 


are spectrally equivalent in the far field. 


3 The impedance Green's function 

We now apply the preceding analysis to the case of the impedance Green's function. The 
impedance Green's function is the solution to the following boundary value problem: 

AGI(x,xo) +k^Gl{x,xo) = S{x-xo), y > 0, 

dG^ (3-1) 

+ ikoc gI — 0, y — 0, 

dn 

where x = {x,y) and xq = (xo,yo)/ with both points located in the upper half-plane. The 
system ( |3.1| l can be re-written as a scattering problem if we express in two parts: 

GI{x,xo) = Gk{x,xo) + u%x), (3.2) 


i.e. the sum of a free-space Helmholtz point-source G^ located at xq and a scattered field 
(see Figure [^. Because of translation invariance in x, we assume that the point source is lo¬ 
cated at xo = {0,h) with h small and positive. 

Let us now denote by O the upper half space y > 0 and by F = 30 the line y = 0. From ( |1.1| | 
and ( |1.2| |, the scattered field w® G G^(n) U G(n) must satisfy 


Aw® + k^u^ = 0 in O, 

3u® „ 

-h ikau = y on F, 

dn 

where x = {x, 0) on F and 

The scattered field w® must also satisfy the radiation condition 


lim i/r 

r^oo 



- fJtw® 


= 0 . 


(3.3) 


(3.4) 
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(3.5) 





The standard approach for evaluating is to write 



^-yJX^-k^y 

VA2 - 


f(A) dA, 


where ^{A) is an unknown function. This represents a solution to the Helmholtz equation, 
enforces the desired radiation condition, and permits the imposition of the boundary condition 
in ( |3.3| > mode by mode. It is straightforward to check that 




a/A^ — + ikoc 

a/— k^ — ikoc 


This is a convenient solution when h is large. When h is small, however, the interval of integra¬ 
tion must clearly be of the order 1 //z, which can be prohibitive. 

While the standard approach to accelerating convergence of the Sommerfeld integral is 
based on contour integration or some variant of the method of images (including complex im¬ 
ages), we seek instead to compute using a combination of the Sommerfeld integral and a 
single layer potential, as we did for the Dirichlet problem above. For this, we let 


= Sro[c^w] +f'io[^w]/ 


(3.6) 


where 

rMo 

SroMW=/ Gk{x,x')cr{x')ds{x'), (3.7) 

J-Mo 

Flo [?] (^) = / (3.8) 

—No y — k^ 

Here, Tq is the finite segment (—Mq, Mq) on the physical interface T and Iq = (—No, No) is a 
finite segment in the Fourier transform domain. We assume crw £ C(ro) and £ C(Io). Note 
that cr(x') = cr(x') with x' = (x',0) on Fo. 

Unlike the Dirichlet problem, we do not have an analytic formula for the single layer density 
cr and cannot determine (Tw by a simple windowing procedure. Instead, we first determine a 
density cr on Fq by letting 

Wl = SfoM 

and enforcing the boundary conditions in ( |3.3| l only on the interval [—Mo, Mo]. Using standard 
jump conditions lUTlITSl , this leads to the local integral equation 

- ^cr + claSroH (3.9) 

on Fo. We will show that cr 6 C(ro) for smooth right-hand-side functions g. After finding cr. 
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we define cr\\i by 


O'w = ^Mo r 


(3.10) 


as above. 


Given cr^y from (3.101, we may substitute it into (|3.6|l and solve for in the Fourier domain; 


fw(A) (-1 + ^ (l + 7^^) (3.11) 


for each A 6 Iq- We discuss the well-posedness of the local integral equation in Section]^ and 
show in Section|^that decays exponentially, in a manner controlled by Mq, the length of the 
window interval Fq, independent of the source location h. 


4 The layered media Green's function 


The scheme described above for impedance boundary conditions can be extended in a straight¬ 
forward manner to the case of layered media. For simplicity, we assume there is a single ma¬ 
terial interface, F, that the wavenumber is fci in the upper layer Qi = {y > 0}, and that the 
wavenumber is k 2 in the lower layer O 2 = {y < 0}. The layered media Green's function is the 
solution to the boundary value problem: 


AGim(x,Xo) + klGim{x, Xo) = S{x - Xo), y > 0, 

AGifn(x,xo) +k2Gim{x,xo) =0, y < 0, 


subject to continuity conditions on T of the form 


[u] = 0, 


du 

dn 


= 0. 


(4.1) 


(4.2) 


We assume that the source lies in the upper half-space at aco = (0, h), radiating at wavenumber 
ki. We then represent the total field in the top and bottom layers as 


Ui{x) = ul{x) + G)tj(x,Xo), for X G Oi, 
W 2 (x) = wKx), for X G 0.2- 


(4.3) 


By analogy with the impedance case, we represent the scattered field in the form 
Ml = [crw] + [yw] + for x G Oi, 

(4.4) 

Ml = s’^] [crw] + d’^I [y w] + F^^ ^r x G O 2 , 

where a and y are unknown charge and dipole densities on Tq. Note that the Sommerfeld 
densities 1 and ^vv 2 are distinct, one invoked for the upper layer and one for the lower 
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layer. As above, we first solve a local integral equation on Fq to obtain functions a and }i. The 
integral equation is derived by enforcing the continuity conditions ( |4.2| ), and classical potential 
theory yields 




(4.5) 


where is the normal derivative of the double layer potential Dp^ on Fq, respectively ITTlITSl . 
Once equation ( |4.5| ) is solved, we let 


cr^Y — Wmo 
flW = ^Mo 


(4.6) 


We then substitute (Tvv and }1\y into the representation ( |4.4| ). Taking the Fourier transform, we 
enforce the continuity conditions (|4.2|l frequency by frequency, and obtain arid ^w,2- 


5 Well-posedness of the integral equation 

Our method relies on the solvability of equations ( |3.9[ | and ( |4.5| ). 

Theorem 5.1. Let gbea Holder continuous function on Tq with exponent a > 0, that is, g G C‘^'“(ro). 
Then there exists a unique solution u G C(ro) to the integral equation 

-^cr + ikaST^,[cr] = g. (5.1) 

Proof. Since Sr^ is a compact operator from C(ro) to C(ro) HTTI . we may obtain the desired 
result by means of the Fredholm alternative and simply show uniqueness for the homogeneous 
case. Thus, assume g = 0 and let 

u(x) = Sro[cr](a;:) forxG]R^\Fo. (5.2) 


A simple calculation shows that v satisfies the boundary value problem 


Av + Cv = 0 


dv+ 

dn 

dv~ 

dn 


+ ikoiv^ — 0 
— ikav~ = 0 


inR^\Fo, 
on To, 
on To, 


(5.3) 
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where 


= lim v(x ± Sn(x)), 


dn 


lim n{x) ■ Vv(x ± Sn(x)), 


(5.4) 


for X 6 Fq. The system < |5.3| l defines an impedance problem on the open arc Fq, which has a 
unique solution Il20l . With zero boundary data, it has only the trivial solution, so that 


dv du+ 
dn dn 


(5.5) 


□ 


Remark 3. It is important to note that in the preceding theorem, Fq is an open interval. The 
density a generally exhibits singularities at the end points ll20ll30H . 

Remark 4. When h > 0, g E C°°(Fo) in equation < |3.9| l and Sr,, is a pseudo-differential operator of 
order minus-one ||25H . This implies that cr = 2{ika SroCr — g) E C‘^'‘*(Fo) with 0 < a < 1. Thus, 
e C~(Fo). 

Remark 5. The proof of existence and uniqueness for the local integral equation in the layered 
media case is analogous and omitted. 


6 Exponential decay of the Sommerfeld integral 


In this section, we outline a proof of the fact that in < |3.11| l is rapidly decaying, independent 
of the source location h. The dielectric case is more involved but the proof follows from the 
same reasoning. 


Theorem 6.1. Let w® be the scattered field induced by a point source at {0,h) satisfying (3.3). Ifzve 
represent w® in the form ( |3.6| ), then ||w('^) | decays superalgebraically, independent ofh. 

Sketch of proof. Rewriting ( |3.11| l slightly, satisfies 

^ ika \ ^ 1 ^ -—:—^ 

-1 + ^^2 _ pi ^ s + 2 ^^ ~ 


— + — ikoi Sfo [cTw]) 


( 6 . 1 ) 


= g-g^ + ’^-ig+\(r-ika Sr, [cr]) + ikx W • [(1 - W)cr] 

— ikoc{\ — W)SrQ[c^w] • 


Each of these terms satisfies the desired superalgebraic decay estimate. The result for g — ^ 
follows the proof of Theorem |2.2[ The second term vanishes since it is the Fourier transform of 
the residual from solving the local integral equation. (In practice, it is zero to discretization er¬ 
ror.) The third term is complicated since the solution of the local integral equation, a, is singular 
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at the endpoints of Fq. However, the frequency content of Sr^ [(1 — W)c 7 '] is controlled near the 
origin, and multiplication by W = Wmq restricts the function to the interval [—3Mo/4,3Mo/4]. 
The last term follows from the spectral equivalence in the far field of Sr^ [c^w] and a mollified 
version of cj-^v. □ 

We state without proof the analogous result for the dielectric case. 

Theorem 6.2. Let and denote the total fields for y > 0 and y < 0, respectively, satisfying the 
continuity conditions ( [4.2| >, with a point source in the upper medium at (0,/z). If we represent w® in the 
form then ||'w,i(A)| and ||^w, 2 (A)| decay superalgebraically, independent of h. 

7 Numerical examples 

In the experiments described below, all of the local integral equations are solved using Nystrom 
discretization on panels, each containing (scaled) 16*-order Legendre nodes. We use gener¬ 
alized Gaussian quadrature for the self-interaction panels, following the approach described 
in ISlIlll. Off-panel evaluation is carried out using Quadrature By Expansion (QBX) 111 [121 
and the resulting linear systems are solved iteratively using GMRES Ii3^ . 

We fix the parameters Mo = 20 and No = 30 for all examples, and as mentioned before, use 
the window function in ( |2.10| ). To accurately evaluate the Sommerfeld integral and avoid the 
square root singularity in the integrand, we deform the integration contour along a hyperbolic 
tangent curve: 

A(t) = t - t e [-30,30]. (7.1) 

This contour is then discretized using 600 uniformly distributed points in t, and the integral is 
evaluated by mean of the trapezoidal rule (which will converge exponentially fast in this case). 

7.1 Impedance Green's function evaluation 

Our first example is simply the evaluation of the impedance Green's function for a source 
located very near to the interface. We set the impedance constant a in equation ( |1.2| ) to be 
a — 1 — O.li. Numerical results are shown in Tablej^and Figure]^ 

Figurej^plots the result when the source is located at (0,10“^) for wavenumber k = 10. One 
can see a sharp spike in the charge density a which causes slow convergence in the Fourier 
domain. However, once we apply the window function to the local integral equation, the 
resulting Sommerfeld density has a rapidly decaying Fourier transform, as shown in Figure [4^ 
By the time A = 30 — ig already less than 10“^®. The density a was discretized along 

To using 160 panels, adaptively refined toward the origin. This required 2560 discretization 
nodes. 
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(a) The real part of the 
impedance Green's function. 
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Figure 4: The impedance Green's function, along with magnitudes of the layer potential 
density a and the Sommerfeld integral density 


Wavenumber k 

1 

1 

10 

10 

20+/ 

20+/ 

Height h 

le-5 

le-8 

le-5 

le-8 

le-5 

le-8 

?(-30) 

7.99e-16 

2.15e-16 

2.06e-ll 

2.06e-ll 

4.83e-ll 

4.83e-ll 

Error 

3.49e-15 

3.28e-15 

2.15e-ll 

2.15e-ll 

1.22e-10 

1.22e-10 


Table 1: Convergence results for the evaluation of the impedance Green's function when 
the source point is close to the interface. 


Table shows the value of the density at f = —30 in ( |7.1| ) and the error in the impedance 
condition du/dn + ikocu at (2.5,0) for a variety of wavenumbers k and heights h. Note that the 
error is independent of the height h and that the density ^ has decayed rapidly 

Remark 6. When the source point is very close to the interface, such ash = 10“^, the maximal 
normal derivative of on T is of the order 0(1 //z), while the potential is of the order 0(log h). 
To avoid catastrophic cancellation, we place an image source at the reflected point across T with 
the same strength as the original source. This cancels the normal derivative of on T. It is 
added back in the final evaluation. We apply the same technique when evaluating the layered 
media Green's function in the next example. This is a finite precision issue, and independent of 
our theory. For a general scattering problem, if the net contribution from a continuous charge 
density on the interface is 0(1), as in Example 3, we do not need to carry out this stabilization. 


7.2 Layered media Green's function evaluation 

In our second numerical example, we evaluate the layered media Green's function. Figure 
shows the real part of the Green's function for ki = 10 in the upper half-space and /C 2 = 20 in the 
lower half-space with source point at (0,10“^). As in the impedance example, there is a spike 
in the dipole density ^ at x = 0. Using only a Sommerfeld integral approach would require a 
prohibitively large interval of integral to obtain convergence, while our hybrid scheme achieves 
rapid convergence in the Fourier domain, as can be seen from the plot of ^vv,i irt Figure]^ More 
detailed data concerning errors in evaluating the Green's functions are provided in Table 
The error is measured as the discrepancy in the potential at (2.5,0) by evaluating the limits 
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(a) The real part of the lay¬ 
ered media Green's function. 



(c) The magnitude of ^w,i- 
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(b) The magnitude of ji. 


Figure 5: The layered media Green's function, the magnitudes of the layer potential density 
fi and the Sommerfeld integral density 


Wavenumber ki 

1 

1 

10+/ 

10+/ 

10 

10 

Wavenumber ^2 

5 

5 

5 

5 

20 

20 

Height h 

le-5 

le-8 

le-5 

le-8 

le-5 

le-8 

ri(-30) 

4.09e-15 

3.67e-15 

1.02e-15 

1.82e-15 

1.63e-ll 

1.63e-ll 

Error 

1.44e-15 

l.Ole-15 

3.88e-14 

3.85e-14 

8.28e-10 

8.28e-10 


Table 2: Convergence results for the evaluation of the layered media Green's function when 
the source point is close to the interface. 


from above and below. The Sommerfeld density decays rapidly, independent of the height /z, 
consistent with our analysis. As with the impedance problem, 160 panels adaptively refined 
toward the origin (2560 points) were required for the discretization of a and fz on Tq. 

7.3 Scattering from inclusions 

We now show the application of our scheme to solving a scattering problem from an object ex¬ 
tremely close to, or partially buried in, a layered media interface. The incident wave is assumed 
to be generated by a point source located at (3,3). 

7.3.1 Close-to-touching object lying in the upper half-space 

We consider a circular object Oq with radius 1 centered at (0 ,1 + e) above the interface. We set 
e — 10“^^, so that the object is only a distance 10“^^ from the interface T. The object is assumed 
to be a perfect conductor or a sound-soft scatterer, both of which require that we impose zero 
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Wavenumber k\ 

1 

1 

10+z 

lO+z 

10 

10 

Wavenumber ^2 

2 

5 

5 

5+/ 

20+5/ 

20 

GMRES iter. 

9 

10 

14 

14 

31 

31 

ri(-3o) 

3.15e-13 

1.08e-13 

3.54e-14 

3.62e-14 

1.78e-8 

4.06e-8 

Error 

9.24e-12 

1.64e-12 

1.32e-10 

2.12e-10 

1.44e-ll 

8.04e-ll 


Table 3: Convergence results for scattering from an object close to the interface of the layered 
media. The closest distance is 10“^^. 


Dirichlet conditions on the object. This yields the following boundary value problem: 


Au + klu — 0, 
Au + = 0, 

u — 0, 
du 


u = 




= 0 , 


in Oi, 
in O 2 , 
on 30o, 

y = 0, 


(7.2) 


along with a suitable decay condition at infinity. 

To solve this scattering problem, we add an additional unknown dipole density yo on the 
boundary of the inclusion Oq. The scattered field w® is then represented by 


1 Sl;W + Df;[)ilYl+F,‘'[fw,ll + D5i^M inn,, 

1 Slj Kl + Dfj [fwI + F,‘= [fw, 2 l in 02 . 

Results for ki = 10 and k 2 — 20 are shown in Figure]^ We only required 120 panels adaptively 
refined toward the origin (1920 points) for the discretization of a and pi on Fq. The discretization 
of }io on OOo required 30 panels. 

From Figure]^ the charge density jio on the boundary of the disk behaves rather benignly. 
However, the density on the interface has a sharp feature near x = 0, as expected. Figure 
shows rapid decay in the Sommerfeld integrand. More detailed errors for various values of the 
wavenumber in the layered media are shown in Table The GMRES iterations are stopped 
when the residual is less than 10“^^. The errors are obtained by solving an artificial scattering 
problem with a known exact solution. For this, the field in each domain is defined by a set of 
free-space sources located in the complement of the domain (see 121 1 for further details). In all 
of our tests, the Sommerfeld integrand decays rapidly. 


7.3.2 Scattering from a partially buried object 

In our final example, a disk of radius 1 crosses the interface between the two layers. The center 
of the disk is (0,0), and we solve the same boundary value problem as in ( |7.2| ). The scattered 
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(a) The real part of the total field scattered in 
the layered media with a circular object lo¬ 
cated right above the interface. 




(b) The real and imaginary parts of the dipole 
density on the boundary of the object 



(d) Real part of the Sommerfeld density ^w,i- 


(c) Real part of the charge density ay^. 

Figure 6: Scattering from an object arbitrarily close to a layered media interface. 



field is represented as: 

si [crw] + Di [fiy,] + Fi [Iwa] + IM 

si [crw] + Di [fiv,] + rf [ew,2] + [^o] 

Note that the density fio on the scatterer is used globally - on both sides of the layered media 
interface. As in the earlier work Il4l l2TH . this has the advantage that the resulting integral 
equation is of the second kind. Results are shown in Figure The densities i and ^w,2 
in the respective Sommerfeld integrals are again rapidly convergent. Along Fq, 180 panels 
(2880 points) adaptively refined toward the intersection of the inclusion and Fq were used to 
discretize cr and /r on Fq. Discretizing pio on OOq required 60 panels. Detailed numerical errors 
are presented in Table The results, again, are consistent with our analysis. 



in Qi, 
in 02 . 


(7.4) 
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Wavenumber k\ 

1 

1 

10+z 

lO+z 

10 

10 

Wavenumber k 2 

2 

5 

5 

5+/ 

20+5/ 

20 

GMRES iter. 

10 

13 

16 

12 

16 

40 

fi(-30) 

7.68e-14 

3.78e-14 

1.88e-14 

1.42e-14 

1.41e-9 

1.32e-6 

Error 

3.47e-12 

7.41e-12 

1.71e-13 

4.07e-13 

5.18e-10 

6.69e-10 


Table 4: Convergence results for scattering from a partially buried object in layered media. 


8 Conclusions 

We have constructed a hybrid approach to acoustic wave scattering in layered media or in half¬ 
spaces with impedance boundary conditions. Our approach retains the advantages of classical 
(physical space) layer potentials in handling close-to-touching interactions and the advantages 
of the Sommerfeld integral in representing smooth interactions along infinite boundaries. By 
solving a local integral equation and subtracting its effect from the original boundary data, we 
have shown that the remaining problem can be solved in the Fourier domain with a rapidly 
convergent integrand. 

We have also shown that the hybrid representation is very convenient when solving scatter¬ 
ing problems with obstacles (including partial buried obstacles). High-order accurate results 
are easily obtained without the direct construction of the Green's function. Instead, our repre¬ 
sentation splits the problem into a free-space scattering problem posed on obstacles and finite 
segments, plus a Sommerfeld correction to enforce the boundary condition along the infinite 
interface. We are currently extending our method to axisymmetric and fully three-dimensional 
problems in both acoustics and electromagnetics. 

Acknowledgements: We would like to thank Alex Barnett, Charlie Epstein, and Tom Hagstrom 
for several useful discussions. 
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Figure 7; Scattering from an object partially buried across a layered media interface. 
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